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Abstract. A new method is presented for mesoscopic simulations of particle 
dispersions in liquid crystal solvents. It allows efficient first-principle simulations 
of the dispersions involving many particles with many-body interactions mediated 
by the solvents. Demonstrations have been performed for the aggregation of 
colloid dispersions in two-dimensional nematic and smectic-C* solvents neglecting 
hydrodynamic effects, which will be taken into account in the near future. 
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1. Introduction 

Dispersions of small particles in host fluids such as colloid suspensions and emulsions are 
of considerable technological importance and appear often in our everyday life as paints, 
foods, and drugs for example. Many kinds of exotic interactions were found between 
particles mediated by the host fluids including screened Coulombic [T|, depletion [I], 
fluctuation induced [2J, and surface induced [Sj forces. A striking example occurs when 
spherical particles are immersed in a liquid-crystal (LC) solvent in the nematic phase 
For a single particle, the orientation of the solvent molecules is distorted due to the 
anchoring of the solvent molecules at the particle surface. Extensive studies have been 
done on this effect, and several characteristic configurations of the nematic field around a 
spherical particle have been identified El HI IH1 E3 E] ■ When the strength of anchoring 
is increased so that normal anchoring is preferred, the solvent changes from quadrupolar 
to dipolar symmetries around the particle. When multiple particles are immersed in the 
solvent, long-range anisotropic interactions are induced between particles due to elastic 
deformations of the nematic field El El E|- The anisotropic interactions can 
have a pronounced effect not only on the local correlations of the particles [H], but 
also on their phase behavior ^SJ El El El El 120] and on their mechanical properties 
[To] . It has been reported also that a similar, but somewhat different situation occurs 
when colloid particles are immersed in smectic-C* films [2U 123 EH] • The purpose of 
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our project is to develop an efficient method suitable for simulating colloid dispersions 
immersed in LC solvents by using smooth interface between colloids and solvents. Effects 
of hydrodynamics is omitted in the present paper though it is very important [2U ESj- 
That will be taken into account in the near future also through the smooth interface 
[27?] . The same method is applicable for charged colloid dispersions where the charge 
density on colloid surface is given by a smooth function [2*Tj . 



2. Mesoscopic Model 

Since analytical approaches for investigating complex materials such as LC colloid 
dispersions considered here are extremely difficult, computer simulations are the 
most promising tool to investigate their static and dynamical properties. In colloid 
dispersions, the host fluid molecules are much smaller and move much faster than the 
dispersed particles. This enables us to use coarse grained mesoscopic variables, which 
should be described by the hydrodynamics, for the host fluids rather than treating 
them as fully microscopic objects [2S1 120] • In the case of charged colloid suspensions 
for example, a mesoscopic method for the first principle simulations can be derived 
by treating the counter ions as a charge density [30] • For the particle dispersions in 
LC solvents considered here, the mesoscopic free energy T of the system can be given 
by functionals of the director n(r), a common direction on which solvent molecules 
are aligned on average with a constraint |n(r)| = 1, for a given particle configuration 
{Ri • ■ ■ Rat}: 

F(n(r);{R 1 ---R N })=F el + F s = J drf el (r) + j> dSf s (S), (1) 

where 

/ e ,(r) = X - [K X (V ■ n) 2 + K 2 (n • V x n) 2 + K 3 {n x (V x n)} 2 ] (2) 

is the Frank free energy jHIj which presents the elastic energy density of the nematic 
solvent at r and the surface free energy 

f s (S) = y[l-(n^) 2 ] (3) 

controls anchoring of the LC solvent at the surface element S of the particle. The 
coefficients Ki, K 2 , and K 3 are the splay, twist, and bend elastic constants, respectively, 
and the single constant approximation (K = K\ = K 2 = K 3 ) reduces Eq.(J2J) to the 
simpler form 

/e*(r) = f [(V • n) 2 + (V x n) 2 ] . (4) 

W is the surface anchoring constant, and v is the unit vector normal to the colloid surface 
[TJ IHJ - The saddle-splay elastic term [H] is not considered here. The integral in the first 
term on the right hand side of Eq.Q runs over the whole solvent volume excluding the 
particles, and that in the second term runs over all solvent-particle interfaces. A simple 
scaling argument tells us T e i oc Ka d ~ 2 and T s oc Wa d ~ l with a and d being the particle 
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radius and the system dimension, respectively, thus the state of the dispersion should be 
controlled by the ratio T s jT e \ oc Wa/K. Although this type of free energy functional is 
sufficient for single particle problems, it is not useful for simulating colloid dispersions 
involving many particles using molecular dynamics (MD) or Brownian type methods 
because the coupling between solvent and the particles is given implicitly by limiting 
the integration space in both T e i and T s . This produces mathematical singularities at 
the interface when one calculates the force, f^ S = ~ &F/dR n , acting on each particle 
mediated by the LC solvents. Calculation of the force is crucial for performing efficient 
simulations of many particle systems. Another serious problem of this type of functional 
is that in order to give correct boundary conditions at the particle-solvent interface, 
one has to use appropriate coordinates for performing grid-based numerical simulations 
rather than the usual Cartesian coordinates. This is generally difficult for particles with 
non-spherical shapes or for systems involving many particles even when each particle has 
a spherical shape. Also this makes the use of the periodic boundary condition difficult, 
which is a fatal situation for simulating bulk materials. 

To overcome these problems, we have modified Eqs.((T]), (J3J), and (jU by using a 
smooth interface between the solvent and the particles so that the coupling is given 



explicitly in the integrand through the interface 
we propose has the form 

^(n(r);{Ri...Rtf}) 
dr 



[32] . The new free energy functional 
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for smectic-C* solvent. The summation convention is used for a, (3, 7 G x,y,z. The 
explicit form of the interfacial profile <p n between dispersed particles and solvents and 
its spatial derivatives V a (f) n are given by 
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Figure 1. The elastic free energy density f e i{%) around a point defect at x = 0. The 
dashed line is from Eq.Q, / e / ~ K/x 2 , and the solid line is from Eq.© or JHJ with 
R c = l, / e , ~iftanh(l/x 2 ). 



respectively with the particle radius a and the interface thickness £. Note that a set of 
Eqs.(l5j)- (fTT|) reduces to Eqs.((T]), Q, and for the limit R c ,(, — > 0. A similar idea for 
using smooth interface was used to treat the hydrodynamic forces acting on particles 
dispersed in simple liquids , recently. For the nematic case, f e i is given by a function 
of a symmetric second-rank tensor q a p = n a np rather than the director n itself to 
take into account the symmetry of the nematic director +n — n automatically. The 
semi-empirical functional form tanh[i?^ • ■ ■] is applied in Eqs.fjHJ) and (JHJ) to avoid the 
mathematical divergence of the elastic free energy density at the defect centers and to 
limit its value to Af e i ~ K/R 2 C and replace its core size with R c as shown in Figure ^ 
Another way to avoid the divergence would be to use the Landau-de Gennes type free 
energy with an order parameter Q a p{r) = Q( r )9a/3( r ), but this requires a prohibitively 
small lattice spacing near the defect points [34J. 

3. Numerical method 

3.1. First-principle Simulation 

The simulation procedure is as follows. 

(i) For a given particle configuration {R 1 ■ ■ ■ Rat}, obtain the interface profile 0„(r) by 
Eq. pOjl . Then we can calculate the stable (or meta-stable) nematic configuration 
n(°) (r) which should satisfy the equilibrium condition 



under the director constraint |n(r)| = 1. One can perform this by numerical 
iterations such as the steepest descent or the conjugate gradient method. Although 
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the time evolutions of n(r) should be determined by the hydrodynamic laws, we 
here assume n(r) follows adiabatically after {Ri • ■ ■ Rat} for the sake of simplicity. 

Once n(°)(r) is obtained, the force acting on each particle mediated by the nematic 
solvents follows directly from the Hellmann-Feynman theorem, 
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(15) 



where Eqs. (ll4j) and (|15j) are for the nematic and smectic-C* solvents, respectively. 
These forms are very convenient because one can compute both d4> n /dR n and 
<9(V a4> n ) I dR n at any time since 4> n is an analytical function of R n . 
Finally, update the particle positions according to appropriate equations of motion 
such as 

gPR*, 



/;i = f PP + fPS + f H + f R 

dt 2 ' 



(16) 



where i pp is the force due to direct particle-particle interactions (hard or soft sphere 
for instance), fjf and fjf are the hydrodynamic and random forces. Repeating the 
steps (i)~(iii) enables us to perform first-principles mesoscopic simulations for the 
dispersions containing many particles without neglecting many-body interactions. 



4. Results 



4-1. Director configurations around a single particle 

We have performed simple demonstrations for two-dimensional (2D) LC colloid 
dispersions to test the performance of our procedure. The demo system has 100 x 100 
lattice sites in a square box with a linear length L = 100. Other physical parameters 
are chosen rather arbitrarily as R c — 1, a — 5, and £ = 2, where the unit of length 
is the lattice spacing I. Since the director configurations in 2D can be expressed by a 
single scalar field 0(r), the tilt angle of the director against the horizontal (a>) direction, 
Eq.((EJ) then reduces to 

dn a {r) bT 

56(r) d9(r) 5n a (r) ' 1 ' 

The boundary condition is fixed at 8(r) = at the edge of the box to avoid rotations of 
the reference frame. 
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I- igurc 2. (a) Director configuration around a single particle immersed in nematic 
solvent with strong anchoring condition Wa/K = 4. The crosses show two —1/2 
charge point defects and the obtained configuration has a quadrupolar character, (b) 
Director configuration around a single particle immersed in smectic-C* solvent with 
Wa/K — 5. The cross shows —1 charge point defect and the obtained configuration 
has a dipolar character. The white disks indicate the colloid particles. Only 9% of the 
total system is shown for display purpose. 



We first calculated the stable director configurations around a single particle 
immersed in 2D nematic and smectic-C* solvents and showed them in Figure El (a) 
and (b), respectively. For the nematic solvent (a), the particle is accompanied by 
two —1/2 charge point defects with a strong anchoring condition Wa/K = 4. The 
distance between the defects and the particle center is about 1.2a, which is similar 
to the analytic value 1.236a [34,, and the director configuration around the single 
particle possess quadrupolar symmetries. This would correspond to the Saturn ring 
configuration observed in three dimensional (3D) systems. Although in principle 
particles can be accompanied by one —1 charge hedgehog defect in 2D as well as in 
3D, such configurations are unstable in the present 2D system since the elastic penalty 
of having m point defects with charge c scales as mKc 2 . This was directly confirmed by 
recent simulations with perfect normal anchoring [33] and also by our simulations. For 
the smectic-C* solvent (b) on the other hand, the particle is accompanied by one — 1 
charge point defect with an anchoring condition Wa/K = 5. The distance between the 
defect and the particle center is about 1.4a which is again similar to the experimental 
value 1.4a±0.1a (23] and the analytic one y2a [2H]- Note that the director configuration 
around the particle possess dipolar symmetries similar to the experiment |21j . 

4-2. aggregation of colloids in LC solvents 

We have simulated the aggregation and ordering process of 30 and 10 colloid particles in 
nematic and smectic-C* solvents, respectively, after the isotropic to nematic or isotropic 
to smectic-C* transition occurred. Here we used the periodic boundary condition and 
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set Wa/K = 4 (nematic) or 5 (smectic C*). Other parameters are the same as in 
the previous single particle case. The simulation was performed starting from a random 
particle configuration which is a typical configuration when the solvent is in the isotropic 
phase (K = 0). We then set K = 1 and calculated f^ s according to the present 
procedure. The particle configurations were updated by numerically solving the steepest 
descent equation, 

c^ = e+fr, as) 

which is obtained by simply substituting d 2 H n /dt 2 — 0, = 0, and = —(dH n /dt 
in Eq.(|16p. ( = 1 is a friction constant and thus the off-diagonal components of the 
hydrodynamic interaction were not considered. Here we obtain 



from the repulsive part of the Lennard- Jones potential, 

N-l N 

E PP = 0.4 £ £ 



2a \ i 2a \ 1 



4 



(20) 



n=l m=n+l 

truncated at the minimum distance |r n — r m | = 2 7 / 6 a, to avoid the particles overlapping 
each other within the colloid radius ~ a. Snapshots from the present simulation 
are shown in Figures H3 and 0] for dispersions in nematic and smectic C* solvemts, 
respectively. In Figure El the particles are forming ordered clusters due to the 
quadrupolar attractive interaction among them. In Figure 0] on the other hand, the 
particles are forming string-like clusters due to the dipolar attractive interaction among 
them. Similar results have been obtained recently by experiments [21]. Note that only 
up to two particle simulations have been done so far [S] and simulations of more than 
two particles would be extremely difficult or almost impossible even for 2D systems by 
means of other methods ever proposed. 



5. Concluding Remarks 

In summary, we have developed a powerful simulation method to investigate colloid 
dispersions interacting via solvents. We proposed a free energy functional which is 
suitable for MD type simulations. The following modifications have been made to the 
original Frank free energy functional. 

(i) The coupling between the nematic solvent and particles at the interfaces is 
introduced explicitly through a smooth interface so that we can analytically 
calculate the force acting on each particle mediated by the host by taking derivatives 
of the free energy according to the particle positions. 

(ii) The value of the free energy density is bounded semi-empirically to avoid a 
mathematical divergence in the defect centers. 




n — / | \ - 

6> 0.5k 7i 

Figure 3. The aggregation and ordering process of colloid particles after the solvent 
exhibit the isotropic (K = 0) to nematic (K — 1) phase transition; (a) t — 0, (b) 
t = 0.4, (c) t = 1.6, (d) t = 6.4, (e) t = 16, and (f) t = 40. Each particle is 
accompanied by two —1/2 charge point defects. Darkness presents the value of q% x - 
Black and white correspond to q^ x — and 0.25, respectively. Those correspond also 
to 9 = 0.257T, 0.757T and 9 = 0, 0.5ir, ir as shown in the gradation map. 
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Figure 4. The aggregation and ordering process of colloid particles after the solvent 
exhibit the isotropic [K = 0) to smectic-C* (K — 1) phase transition; (a) t = 0, (b) 
t = 0.25, (c) t = 1, (d) t = 4, (e) t = 10, (f) t = 25. Each particle is accompanied 
by one —1 charge point defect. Darkness presents the value of q% x - Black and white 
correspond to q% x = and 0.25, respectively. Those correspond also to 8 = 0.257T, 0.757T 
and 9 — 0, 0.57T, tt as shown in the gradation map. 
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We have performed demonstrations for 2D dispersions and confirmed that the method 
works quite well for systems which contain defects. Applications of this method 
to 3D systems should have no theoretical difficulties, but require somewhat heavier 
computation. This should allow the simulation of the chaining of the particles caused 
by the possible dipolar symmetry of the nematic configurations around a single particle. 
Although we have shown only simple demonstrations of the method by performing 
simulations of the 2D system in this letter, simulations with physically more interesting 
situations such as systems with non-circular particles, asymmetric particle pairs with 
different particle size, or particles with non-normal anchoring as well as more realistic 
simulations in 3D systems are possible. Efforts for taking into account the hydrodynamic 
effects are now underway by using smooth interface instead of imposing hydrodynamic 
boundary conditions at the colloid interface. 
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Appendix A. numerical implementation for 2D nematic 

Here, the superscripts i,j (— 1, 2, • • • , L) denote positions on a 2D lattice (L x L). 
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Appendix B. numerical implementation for 2D smectic-C* 
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